In [12]:

    

# Let's check our version
!cs2cs


    

    




Rel. 4.7.1, 23 September 2009
usage: cs2cs [ -eEfIlrstvwW [args] ] [ +opts[=arg] ]
                   [+to [+opts[=arg] [ files ]

In [13]:

    

# We would expect this to work, but it doesn't
# The towgs84 is used for transformation
! echo "49554 215020 0" | cs2cs -v +init=epsg:31370 +to +init=epsg:4326  -f %.6f


    

    




# ---- From Coordinate System ----
#Lambert Conformal Conic
#	Conic, Sph&Ell
#	lat_1= and lat_2= or lat_0
# +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339
# +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438
# +ellps=intl
# +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m
# +no_defs
# ---- To Coordinate System ----
#Lat/long (Geodetic alias)
#	
# +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
# +towgs84=0,0,0
#--- following specified but NOT used
# +ellps=WGS84
2.928019	51.235804 347.868509

In [14]:

    

# We get a different result if we use invproj, what's going on?
! echo "49554 215020" | invproj -v +init=epsg:31370 -f %.6f


    

    




#Lambert Conformal Conic
#	Conic, Sph&Ell
#	lat_1= and lat_2= or lat_0
# +init=epsg:31370 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339
# +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438
# +ellps=intl
# +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m
# +no_defs
2.929249	51.236883

In [15]:

    

# We get the correct result if we fix the last parameter of the towgs84 parameter
# http://osgeo-org.1560.x6.nabble.com/proj4-epsg-31370-possible-wrong-parameters-td3841364.html
!echo "49554 215020" | cs2cs -v \
    +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 \
    +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl \
    +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m \
    +no_defs -to  +proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs -f %6f


    

    




o ---- From Coordinate System ----
#Lambert Conformal Conic
#	Conic, Sph&Ell
#	lat_1= and lat_2= or lat_0
# +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
# +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747
# +units=m +no_defs +datum=WGS84
#--- following specified but NOT used
# +proj=longlat +ellps=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
o ---- To Coordinate System ----
#Lat/long (Geodetic alias)
#	
# +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
2.930480	51.236358 41.428643

In [16]:

    

# This does not work for invproj because proj and invproj ignore the datum shift
# The towgs84 is ignored
!echo "49554 215020" | invproj -v \
 +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 \
 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl \
    +towgs84=-1,1,1,1,1,1,1 +units=m \
 +no_defs \
 -f %6f


    

    




#Lambert Conformal Conic
#	Conic, Sph&Ell
#	lat_1= and lat_2= or lat_0
# +proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
# +towgs84=-1,1,1,1,1,1,1 +units=m +no_defs
2.929249	51.236883

In [17]:

    

# There is another datums shift that gives the same result on spatialreference.org
# http://spatialreference.org/ref/sr-org/epsg31370-correct-datum-shift/
!echo "49554 215020" | cs2cs -v \
+proj=lcc +lat_1=51.16666723333334 +lat_2=49.83333389999999 \
+lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 \
+ellps=intl +towgs84=-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0 +units=m +no_defs \
-to +datum=WGS84 +ellps=WGS84 -f %6f


    

    




o ---- From Coordinate System ----
#Lambert Conformal Conic
#	Conic, Sph&Ell
#	lat_1= and lat_2= or lat_0
# +proj=lcc +lat_1=51.16666723333334 +lat_2=49.83333389999999 +lat_0=90
# +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
# +towgs84=-99.1,53.3,-112.5,0.419,-0.83,1.885,-1.0 +units=m +no_defs
# +datum=WGS84
#--- following specified but NOT used
# +ellps=WGS84 +ellps=WGS84 +towgs84=0,0,0
o ---- To Coordinate System ----
#Lat/long (Geodetic alias)
#	
# +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
2.930478	51.236358 41.262401

In [18]:

    

%%file test.vrt
<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>test.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:31370</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="x" y="y"/>
    </OGRVRTLayer>
</OGRVRTDataSource>


    

    




Overwriting test.vrt

In [19]:

    

! echo "x,y\n49554,215020"  >test.csv 
jsontxt = ! ogr2ogr -f GeoJSON -t_srs EPSG:4326 /dev/stdout test.vrt
import json
#json.loads(jsontxt)
print("\n".join(jsontxt))
json.loads("".join(jsontxt))['features'][0]['geometry']['coordinates']


    

    




{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "x": "49554", "y": "215020" }, "geometry": { "type": "Point", "coordinates": [ 2.930479615690256, 51.236357952698896 ] } }

]
}






    
Out[19]:





[2.930479615690256, 51.236357952698896]

In [20]:

    

# Or you can use a python script to talk to the osr library, which as the correct coordinates.
import osgeo.osr
src = osgeo.osr.SpatialReference()
dst = osgeo.osr.SpatialReference()
src.ImportFromEPSG(31370)
print src.GetTOWGS84()
dst.ImportFromEPSG(4326)
transform = osgeo.osr.CoordinateTransformation(src, dst)
lon, lat, z = transform.TransformPoint(49554, 215020, 0)
print(lon, lat, z)


    

    




(-106.869, 52.2978, -103.724, 0.3366, -0.457, 1.8422, -1.2747)
(2.930479615690256, 51.236357952698896, 41.428643393330276)

In [20]: